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ABSTRACT 


A  new  two-dimensional  data-adaptive  algorithm  utilizing  the  iterative  Toeplitz 


\ 


approximation  method  is  presen^d.  This  algorithm  provides  a  robust  and  efficient 
means  for  accurate  estimation  of -2-^  autoregressive  parameters  representing  spatially 
variant  data  arrays.  Its  convergence  performance  is  comparable  to  that  of  the  2-D 
Recxirsive  Least  Squares  (RLS)  algorithm  but  is  computationally  more  eflBcient.  The 
savings  in  computation  is  realized  by  reducing  the  size  of  the  matrix  to  be  inverted 
when  solving  the  AR  model  normal  equations.  The  ability  of  the  algorithm  to  accu¬ 
rately  estimate  the  model  parameters  using  very  small  data  sets  and  various  window¬ 
ing  schemes  are  evaluated.  Spectral  estimates  are  compared  for  quarter-plane  (QP), 
nonsymmetric  half-plane  (NSHP)  and  combined-quadrant  (CQ)  regions  of  support. 
Additionally,  the  algorithm  is  tested  in  noise  cancellation  and  line  enhancement  ap¬ 
plications  using  image  arrays.  This  algorithm  may  be  implemented  for  data-adaptive 
image  processing  or  coding  and  for  applications  requiring  2-D  spectral  estimation. 
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I.  INTRODUCTION 


A.  PARAMETER  ESTIMATION  OF  DATA-ADAPTIVE  ARRAYS 

Two-dimensionaJ  (2-D)  autoregressive  (AR)  modeling  provides  a  means  for  the 
spectral  estimation  of  signals  received  by  an  array  of  spatially  distributed  sensors. 
Additionally,  the  2-D  AR  model  parameters  may  be  used  to  represent  the  original 
signal  in  a  compact  manner  resulting  in  compression  of  the  data  to  be  processed  or 
transmitted.  This  is  particularly  useful  in  many  image  processing  applications.  The 
traditional  least  squares  or  Wiener  filtering  based  2-D  AR  parameter  estimation  [Ref. 
1]  provide  accurate  spectral  estimates  but  require  the  the  inversion  of  a  relatively  large 
autocorrelation  matrix.  This  is  undesirable  since  matrix  inversion  is  computationally 
expensive.  Iterative  methods,  using  a  Toeplitz  approximation  algorithm,  have  been 
suggested  [Ref.  2,  3]  as  an  alternative  to  the  direct  Wiener  filter  solution.  The  itera¬ 
tive  Toeplitz  approximation  method  has  proven  to  be  an  effective  means  to  estimate 
the  parameters  of  fixed  data  arrays  without  having  to  invert  the  correlation  matrix. 
This  method  is  particularly  useful  when  data  arrays  of  limited  size  must  be  processed. 

The  results  obtained  using  iterative  methods  to  obtain  the  AR  p2U’ameters  of 
fixed  data  arrays  suggest  that  these  methods  may  be  used  to  some  adveintage  when 
estimating  parameters  adaptively.  Adaptive  parameter  estimation  is  necessary  in 
many  digital  signal  processing  applications  where  the  data  is  non-stationary.  In  these 
applications,  real-time  or  near  real-time  processing  is  often  desirable.  This  requires 
that  <iny  algorithm  used  must  be  computationally  efficient  and  be  capable  of  providing 
accurate  estimates  obtained  from  a  minimum  of  data.  This  thesis  addresses  the 
implementation  of  the  iterative  Toeplitz  approximation  method  for  2-D  parameter 
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estimation  of  uon-stationary  data.  It  is  shown  that  this  algorithm  cam  be  efficiently 
implemented  to  obtain  accurate  estimates  from  very  smzdl  data  arrays. 

B.  EXPLANATION  OF  NOTATION 

All  rectangular  matrices  are  denoted  by  a  boldface,  captial  letter,  e.g.,  R.  Col¬ 
umn  vectors  are  designated  by  a  boldface  lowercase  letter,  for  example:  x.  Occasion¬ 
ally,  a  vector  created  temporarily  for  the  development  of  a  mathematical  expression 
will  be  represented  by  a  lowercase,  boldface,  greek  letter.  An  example  of  this  would 
be  the  vector  a .  A  special  operator  matrix  will  be  represented  by  a  calligraphic  cap>- 
ital  letter,  such  as  J-.  Scalar  values  are  generally  represented  non-boldface  lowercase 
arabic  or  greek  letters,  e.g.,  a  or  A.  Captital  letters  such  as  A  or  P  are  used  to  rep¬ 
resent  a  variety  of  values  or  expressions.  These  include  anything  from  a  polynomial 
to  the  dimension  of  a  matrix.  The  symbol  T  *s  reserved  to  indicate  the  trmspose  of 
a  vector  or  matrix. 

C.  THESIS  OUTLINE 

The  following  describes  the  organization  of  the  remainder  this  thesis.  Chapter 
II  serves  as  an  introduction  to  the  problem  of  2-D  AR  modeling  and  parameter  esti¬ 
mation.  Particular  attention  is  given  to  the  formulation  of  the  normal  equations  for 
quarter-plane  (QP)  and  nonsymmetric  half-plane  (NSHP)  regions  of  support.  This 
chapter  concludes  with  a  description  of  the  least  squares  algorithm  and  direct  in¬ 
version  method  of  solving  the  normal  equations.  This  serves  as  the  foundation  for 
Chapter  III  which  extends  the  concept  of  solving  normal  equations  using  the  itera¬ 
tive  Toeplitz  approximation  method  to  derive  parameter  estimates  referred  to  as  the 
fixed  data  method.  Chapter  IV  proceeds  with  the  development  of  the  data-adaptive 
algorithm  using  Toeplitz  approximation.  Two-dimensional  AR  spectral  estimation 
is  discussed  in  Chapter  V,  along  with  a  comparison  of  results  using  both  fixed  and 
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adaptive  iterative  methods  for  spectral  estimation.  Additionally,  a  discussion  of  the 
use  of  the  data-adaptive  algorithm  for  noise  cancellation  and  line  enhancement  of 
image  arrays  is  provided  in  Chapter  VI,  with  experimental  results.  Conclusions  and 
recommendations  for  future  work  are  found  in  Chapter  VII. 
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II.  2-D  AUTOREGRESSIVE  MODELING 


A.  OVERVIEW 

This  procedure  assumes  a  stationary  (homogeneous)  random  process  x[ni,n2] 
that  is  the  response  of  an  AR  model  excited  by  a  white  noise  input  tn[ni,n2]  having 
a  variance  The  AR  model  as  shown  in  Figure  2.1  may  be  described  as  an  all  pole 


w{ni,n2) 


x(ni,n2) 


Figure  2.1:  AR  model  excited  by  white  noise, 
filter  with  the  transfer  function 

“  1  +  E  AM)  o(i., 

where  A  is  the  region  of  support  over  which  the  parameters  a{ki,k2)  are  non-zero. 
The  difference  equation  for  the  system  that  generates  x[ni,n2]  can  be  written  as 

x[ni,n2]  =  -  *"  *->”2  -  j]  +  w[ni,n2].  (2.2) 

(»j)  €A 

A  linear  set  of  equations  for  the  filter  coefficients  a(z,  j)  may  be  formed  by  multiplying 
both  sides  of  (2.2)  by  z[ni  —  /i,n2  —  /2]  and  computing  the  the  statistical  expectation 
of  the  resulting  expression  [Ref.  4].  This  leads  to  the  following  expression  called  the 
normal  equation 


y  =  -  E  E  <■(;.»  A|(,  -i,i2-  j] , 

(«.i) 


4 


for  li,l2  >  0.  The  coefficients  a{i,j)  can  be  derived  from  this  normal  equation.  It 
should  be  noted  that  the  structure  of  the  normal  equation  depends  on  the  region  of 
support  A.  A  region  of  support  may  have  any  shape,  but  th€^  most  commonly  used 
are  the  quarter-plane  and  nonsymmetric  half-plane. 

1.  Quarter-Plane  Support 

The  most  straightforward  region  of  support  is  the  quarter-plane.  A  region 
A  is  considered  to  have  quarter-plane  (QP)  support  when  a{i,j)  has  non-zero  values 
in  one  quadrant  only  as  shown  in  Figure  2.2.  For  this  case  the  normal  equation  has 


Figure  2.2;  The  quarter-plane  region  of  support. 


the  form 


R.\hM  =  E  (i,i)  (0,0) 

i=0  j=0 


(2.4) 


where  /i  =  1, 2, . . . ,  Pi  — 1  and  =  1,2, . . . ,  Pj— 1  with  Pi  and  P2  being  the  dimensions 
of  A.  If  it  is  assumed  that  a(0,0)  =  1,  we  may  express  (2.4)  in  block  matrix  foriln  as 


r  Ro 

R-i 

R-2 

R 

Ri 

Ro 

R-1 

' "  ’  R 

R2 

Ri 

Ro 

R 

.Rpi-i 

Rpi-2 

RPi-3 

.  .  . 

-Pi+i 

-Pi +2 

-Pi +3 

Ro  J 


ao 
ai 
a2 

Lap,-i  J 


£:(o) 

0 

0 


(2.5) 
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Each  block  Ra/  of  this  matrix  is  given  by 

Rat  —  I^-Af  ~ 


R{M,0) 

R{M,1) 


RiM,-l) 

R{M,0) 


R{M,-P2  +  1) 
R(M,  -Pi  +  2) 


(2.6) 


R{M,P2-\)  R(M,P2-2)  •••  R(M,0) 

while  the  coefficients  are  hm  =  (a(M,0)  a(M,  1)  a(M,2) . . .  a(Af,  P2  —  1)]^  and  the 
error  variance  vector  is  =  [af  0 . . .  0]^.  The  blocks  Ra/  along  the  diagonals  of  the 
autocorrelation  matrix  R  are  equal  and  the  diagonal  elements  of  Ra/  are  equal,  thus 
R  is  Toeplitz  block- Toeplitz. 

2.  Nonsymmetric  Half-Plane  Support 

A  region  A  with  non-zero  a(t,j)  as  shown  in  Figure  2.3  is  considered  to 
have  nonsymmetric  half-plane  (NSHP)  support.  The  normal  equation  for  this  case 


must  be  modified  to  account  for  non-zero  a(0,j)  in  the  fourth  quadrant.  This  may 
be  written  as 

fo+Pj-l  Pj-lPa-l 

J2I  =  <*(0,j)Rxlli,l2  -i]  +  53  53  -  i,l2-jJ  ,  (2.7) 

i=l  »=0  jsrO 
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where  (t,  j)  ^0,  /i  =  0, 1, 2, . . . ,  Pi  —  1  and 


_  f  0?  2, . . . , 

\  i>2,L2  +  1, 


L2  +  P2  —  1> 


for  /i  =  0 


L2  +  2, . . . ,  Ij2  d"  Pi  —  for  ^1^0 


The  dimensions  of  A  axe  given  by  Pi  and  P2  and  L2  is  a  negative  number  as  defined 
in  Figure  2.3.  As  in  the  quarter-plane  case,  if  we  let  a(0,0)  =  1  we  may  write  (2.7) 
in  block  matrix  form  as 


Rk)  R-1  R-2 

Ri  Ro  R-i 

R2  Ri  Ro 

.Rp,_i  Rpi-2  Rpi-3 


R-Pi-j-i  So 
R-Pi+2 

R-Pj+3  ^2 
Ro  L®Pi-l, 


-£(0)- 

0 

=  0 


The  blocks  Ro,  Ra/  and  Km  each  have  different  forms.  The  diagonal  block 


Ro  has  the  form 


Ro  = 


P(0,0) 

R(o,i) 


P(0,-1) 

P(0,0) 


P(0,/2  +  P2-l)  P(0,X2  +  P2-2) 


P(0,  —Li  ~  P2  -f  1) 
P(0,  —Li  —  P2  +  2) 

P(0,0) 


the  blocks  Km  =  RIm  in  the  first  row  and  column  maj'  be  written  as 


R{M,-L2) 
R{M,-L2  +  1) 


R{M,  -L2  - 1) 
R{M,  -L2) 


Lp(M,-I2  +  P2-1)  P(A/,-T2  +  P2-2)  ••• 

and  the  remaining  blocks  Rm  are  given  by  (2.6). 

The  model  parameter  vectors  are  given  by 


P(M,-L2-P2  +  1)' 
P(M,-L2-P2  +  2) 

P(M,0) 


,  (2.10) 


ao  =  [a(0, 0)  a(0, 1)  a(0, 2) . . .  a(0, 12  +  P2  -  1)]’ 


=  [a(M,  L2)  a{M,  L2  +  I)  a(M,  L2  +  2)...  a(M,  L2  +  P2-  1)]^ 
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for  Af  =  1, 2, . . . ,  Pi  —  1.  The  error  vaxiance  vector  6^°^  =  [cr*  0 . . .  0]^.  Except  for  the 
upper  diagonal  block  (2.9)  and  top  amd  left  most  blocks  (2.10),  the  NSHP  autocorre¬ 
lation  matrix  is  block-toeplitz  with  each  block  being  toeplitz  as  well.  Quarter-Plane 
support  can  be  considered  a  special  case  of  NSHP  support  with  L2  =  0. 

B.  SOLVING  THE  NORMAL  EQUATIONS  BY  DIRECT  INVERSION 
The  normal  equations  also  arise  in  the  2-D  linear  prediction  problem.  When 
solving  the  normal  equations,  the  objective  is  to  find  the  parameters  which  mini¬ 
mize  the  prediction  error  [Ref.  5].  It  is  frequently  more  convenient  to  describe  the 
formulation  of  the  normal  equations  in  terms  of  linear  prediction. 

In  the  2-D  least  squares  problem,  the  error  between  the  random  process  r  [ni, 
and  its  estimate  x[ni,n2]  is  given  by 

e[ni,n2]  =  x[ni,n2]  -  x[ni,n2]  (2.11) 

or  in  expanded  form  as 

e[ni,n2]  =  x[ni,n2]  +  532I®(*>iM”i-*»”2-j],  (t,  j)  (0,0).  (2.12) 

(ij)  eA 

The  objective  of  the  least  squares  method  is  to  minimize  the  squared  errors  from  a 
particular  set  of  these  terms  [Ref.  4].  If  we  let  a(0,0)  =  1,  we  can  express  (2.12)  as 

e[ni, 712]  =  ”  jJ  (2.13) 

(»j)  Ci* 

which  is  compactly  represented  in  vector  notation  as  e  =  Xa.  The  error  e[ni,n2] 
is  computed  for  each  position  of  the  filter,  then  normal  equations  are  formed  by 
multiplying  both  sides  of  (2.13)  by  X^  and  applying  the  orthogonality  principle  [Ref. 
5].  This  results  in 

Ra  =  S  (2.14) 
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where  R  =  X^X  and  the  error  vaxiance  vector  £  =  X^e  =  0, . . . ,  0]^  with 

=  [<Tj ,  0, . . . ,  0]^.  The  matrix  X  may  be  defined  as 


■  1 

r  1  1  11 

Xo  X' 

i 

= 

Xo  Xi  •••  XPj_i 

1  1  1 

L  1  J 

L  1  1  1  J 

(2.15) 


and  the  model  paxameters  may  be  given  by  a  =  [1  a']^.  The  parameter  estimates 

are  obtained  from  (2.14)  by  muliplying  both  sides  of  the  equation  by  the  inverse  of 
R  which  may  be  written  in  expanded  form  as 


a'  =  -(X'^X')~'X'^Xo  . 


(2.16) 


This  is  referred  to  as  the  direct  inversion  method. 

The  preceding  discussion  assumes  that  R  is  Toeplitz  block-Toeplitz.  This  true 
for  the  case  where  the  autocorrelation  method  [Ref.  5]  is  used  to  form  the  correlation 
matrix  or  when  the  theoretical  correlation  is  known.  In  practicality,  however,  param¬ 
eter  estimates  must  be  obtained  from  relatively  small  sets  of  finite  data.  In  this  case, 
it  is  often  desirable  to  compute  R  using  the  covariance  method.  This  results  in  a 
block  autocorrelation  matrix  R  with  non-Toeplitz  blocks  Ra/  [Ref.  2].  The  rows  of 
the  matrix  X  are  the  elements  of  the  random  process  covered  by  the  filter  mask  for 
each  point  being  estimated  as  the  mask  is  moved  over  the  data.  In  the  covariance 
method  the  filter  mask  is  not  moved  past  the  boundaries  of  the  region  of  support. 
This  means  that  when  QP  support  is  used,  X  will  be  of  dimension  P1P2  x  P1P2  where 
Pi  and  P2  are  the  dimensions  of  the  data  array  accessible  to  the  filter  mask.  For 
example,  using  a  3  x  3  (nine  element)  filter  meisk  to  form  the  normal  equations  of  an 
8x8  data  array  would  result  in  an  X  matrix  with  dimension  36  x  9.  It  should  be 
noted  that  Pj  =  P2  =  6,  since  two  rows  and  two  columns  of  data  cannot  be  estimated 
by  the  filter. 

The  preceding  discussion  provides  the  basis  for  the  subsequent  work  in  this 
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thesis.  The  next  chapter  will  addresses  this  same  problem  using  the  iterative  Toeplitz 
approximation  method. 
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III.  ITERATIVE  TOEPLITZ 
APPROXIMATION  ALGORITHM 


The  major  drawback  of  the  least  squares  method  described  in  the  previous 
chapter  comes  from  the  requirement  to  invert  a  relatively  large  autocorrelation  ma¬ 
trix.  This  is  undesireable  since  matrix  inversion  is  computationally  expensive.  In 
this  chapter  we  present  an  iterative  method  for  solving  the  normal  equations  which 
does  not  involve  the  direct  inversion  of  the  correlation  matrix  to  obtain  the  model 
parameters.  This  method  is  then  extended  to  the  case  where  the  covariance  method 
is  used  to  form  the  correlation  matrix  and  the  normal  equations  are  solved  iteratively 
using  Toeplitz  approximation. 

A.  THE  ITERATIVE  SOLUTION  OF  THE  NORMAL  EQUATIONS 
An  alternative  to  direct  inversion  is  to  take  advantage  of  the  near  Toeplitz- 
block- Toeplitz  structure  of  the  autocorrelation  matrix  and  successively  partition  the 
normal  equation  [Ref.  2].  This  partitioning  permits  the  estimation  of  parameters 
while  requiring  the  inversion  of  a  matrix  that  is  significantly  smaller  than  the  original 
autocorrelation  matrix.  The  following  discussion  develops  the  iterative  solution  for 
QP  and  NSHP  regions  of  support. 

1.  Quarter-Plane  Support 

In  order  to  iteratively  solve  the  QP  normal  equations  we  must  begin  by 
dividing  both  sides  of  (2.5)  by  <7^.  This  results  in  a  modified  normal  equation  given 
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by 


Ro  R_i  R_2  •  •  R-p,+i 

■g{0)- 

Ri  Ro  R— 1  R- Pi+2 

< 

0 

R2  R-l  Ro  R-Fi+3 

= 

0 

•  1  •  •  • 

.Rp,_i  Rpj_2  Rp,-3  •••  Ro  . 

0 

where  f  =  (1 0  0, . . . ,  0]^  and  aj'  =  ^aj  for  *  =  1, 2, . . . ,  Pi  —  1.  We  now  partition 
the  normal  equation  (3.1)  as  follows 


fGi 

hi] 

[  1  _  1 

^  ll 

(Q  o’i 

W 

RoJ 

Up-i\  I 

0  J  ’ 

where 

r  Ro 

R_ 

1  R_2 

•••  R_p,+2" 

Ri 

Ro 

R-1 

•••  R-P,+3 

Gi  = 

R2 

Ri 

Ro 

R-Pi+4 

,  (3.3) 

.  Rpi  -2 

Rp,. 

-3  Rpi-4 

*  Rq 

=  [Rpi 

_i  Rpi_2  Rpj-j 

Ri]  » 

(3.4) 

ai 

=  (a«'  a;^ 

•f  •• 

•  api-2J  > 

(3.5) 

and 

IT  1  = 

,[^{0)T  qT 

o^f  . 

(3.6) 

We  may  now  form  a 

set  of 

coupled  iterative  equations  by  defining  an  op- 

erator  Pi 

which  is  given 

as 

^>  =  i 

0 

0  R^' 

y 

(3.7) 

and  using 

this  operator  to  premultiply  both  sides 

of  (3.2)  to  yield 

01 

=  Gr' 

[t  1  -  hi  a 

"  1 

Pi-iJ  > 

(3.8) 

and 

4,-1 

=  -R^'hfa,  . 

(3.9) 
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Equations  (3.8)  and  (3.9)  suggest  an  iterative  solution  and  can  be  written  as 

al‘>  =  Gr’h.  -I>,  (3.10) 


and 


*P,-1  —  “l  “l 


(3.11) 


which  provide  a  means  to  solve  the  normal  equation  (3.1)  iteratively. 

It  can  be  noted  that  the  submatrix  Gi  is  nearly  the  same  dimension  as 
the  original  correlation  matrix.  This  means  that  the  inversion  of  Gi  provides  little 
advantage,  computationally,  over  the  direct  inversion  method.  It  is  clear,  however, 
that  further  partitioning  of  Gi  will  decrease  the  required  complexity.  Partitioning  of 
Gi  yields 


where 


'G2 

h2l 

'  1  _  r-y  2I 

[hj 

ReJ 

la^.J  0  J 

r  Ro 

R.1 

R-2 

•••  R_Pj+3‘ 

Ri 

Ro 

R-1 

•••  R-P,+4 

II 

0 

R2 

Ri 

Ro 

•  R-P,+5 

.RPi-3 

RPl-4 

Rpl-5 

Ro  . 

hj—  [Rp,_2  R-P)-3  RPi-4 


Ri  ]  > 


(3.12) 


(3.13) 


(3.14) 


2md 


^  r  ^"T 


_  _  -If 

q2  =  laj) 


•••  “Pi-sJ  > 


(3.15) 


T,  =  [£('»’■  0^  ...  Ol'f  .  (3.16) 

As  before,  both  sides  of  (3.12)  are  premultiplied  by  an  operator  defined  by 

Gr' 


:F2  = 


2^  0 
0  ’ 


(3.17) 
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which  results  in 

02  =  Gj  2  “  ^2  ®Pj_2]  1  (3.18) 

and 

a'A-2  =  -  (3.19) 

Equations  (3.18)  and  (3.19)  may  be  solved  iteratively  in  the  same  m2umer  as  (3.8) 
and  (3.9).  This  solution  is  expressed  as 

a  5*^  =  1  -  ha  a^i~2^^]  »  (3-20) 


and 

4S-2  =  .  (3.21) 

This  partitioning  process  may  be  continued  until  Gp,_i  =  Ro.  This  will 
result  in  a  partitioned  normal  equation  given  by 


where  =  Ri,  ap,_i  =  <  and-y 

Premultiplication  of  both  sides  of  (3.22)  by  zm  operator  Tp^-i  will  result 
in  the  iterative  solution 


a?‘>  =  Rj'fW  -  [R-,af +  . . .  +  ,  (3.23) 

and 

=  -R^»[R,a.?*-')  +  R_ia?*-*>  +  . . .  +  R_p,+2apSiV)] .  (3.24) 

For  this  case  the  operator  matrix  is  given  by 

^''-=[T  vl-  p-25) 
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At  this  point  in  our  development,  we  can  combine  the  preceding  iterative 
solutions  from  the  successive  levels  of  partitioning  into  a  set  of  compact  iterative 
summations 

=  ai'l”’  -  V  ,  (3.26) 

»=1 

af  =  -H.-  s'  R,-.a,"'‘-> ,  (i  #  i) ,  (3.27) 

t=0 

where  and  j  =  1,2, . . . ,  Pi  —  1.  The  index  of 

iteration  is  k  and  aj'  are  the  Pj  x  1  parameter  vectors.  The  solution  of  (3.26)-(3.27) 
requires  OiP^P^K)  multiplications,  where  K  represents  the  number  of  iterations 
to  converge  to  the  true  psjrameter  values.  An  additionzd  P1P2  multiplications  are 
required  to  compute  the  initial  parameter  estimates  a”^°^  for  j  =  1,2, . . .  ,P2  —  1. 
This  results  in  total  multplications  of  OIP^P^K  +  PiP^)-  Depending  on  the  number 
of  iterations  required  to  converge  to  the  true  parameters,  this  total  compares  with  the 
algorithms  of  Wax  and  Kaiaith  [Ref.  6]  and  Akaike  [Ref.  7]  which  require  C7(Pf  P2 ) 
operations.  When  large  arrays  are  used  this  represents  a  noticeable  savings  over  direct 
inversion  which  requires  0{PiI^)  multiplications. 

2.  Nonsymmetric  Half-Plane  Support 

The  derivation  of  the  iterative  NSHP  solution  follows  a  procedure  of  suc- 
cesive  partitioning  similar  to  that  of  QP  support.  However  some  differences  exist  due 
to  the  tisymmetry  of  the  autocorrelation  matrix  consisting  of  three  different  types 
of  block  matrices.  This  asymmetry  results  in  a  somewhat  modified  set  of  recursive 
summations  given  by 

-(*)  ^  _  ^-1  ,3  28) 

>sl 

a"'*'  =  -R,-'R,--ai"‘-‘'  -  R^'  Z  (.'  #  ;),  (3.29) 
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where  and  j  =  1,2, . . . ,  Pi  —  1.  The  index 

of  iteration  is  k  and  aj'  are  the  Pj  x  1  parameter  vectors.  As  with  QP  support.  The 
solution  of  (3.28)-(3.29)  requires  OiPiP^K)  multiplications,  where  K  represents  the 
number  of  iterations  to  converge  to  the  true  psu'suneter  values. 

B.  THE  ITERATIVE  SOLUTION  USING  TOEPLITZ  APPROXIMA¬ 
TION 

The  discussion  of  the  previous  section  assumes  an  autocorrelation  matrix  with 
Toeplitz-block-Toeplitz  structure.  In  most  cases  the  autocorrelation  matrix  must  be 
formed  from  a  limited,  and  possibly  very  small,  amount  of  sampled  data.  In  this 
situation,  it  is  best  to  use  the  covariance  method  to  estimate  the  correlation  matrix 
to  reduce  the  bias  in  the  estimate.  However,  since  this  method  does  not  have  Toeplitz- 
block'Toeplitz  structure,  it  is  necessary  to  modify  the  iterative  algorithm  described 
above. 

1.  Forward  Iteration 

For  first  quadrant  QP  support,  the  covariance  estimate  R  has  the  form 


Ro,o 

Ro.i 

Ro,2 

Ro.Pi-1 

Rin 

Ri.i 

Ri,2 

Ri,/»,_i 

^■2,0 

1^2,1 

R-2,2 

R-2,Pi-1 

(3.30) 

B-Pi-i.o 

R-P,-1.2 

*  •  1 
•••  Rp,-i,p,_i. 

A  proven  procedure  [Ref.  8]  that  can  be  used  to  form  the  Toeplitz  approx¬ 
imation  is  to  first  average  the  diagonal  blocks  Rjj,  *  =  j  and  then  form  a  Toeplitz 
T  matrix  from  this  averaged  block  Ra„,  by  averaging  its  diagonal  elements.  The 
diagonal  elements  of  T  can  be  obtained  as  follows 

J  fli-i-l 

— 7  E  (3.31) 

where  »  =  0, 1, . . . ,  Pj  -  1. 
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We  may  now  proceed  to  succesively  partition  the  normal  equations  in  the 
manner  outlined  in  the  previous  section.  The  first  partitioning  leads  to  the  iterative 
equations 

=  Gf'l-r  .  -  h,  ,  (3.32) 


and 


>"(*)  —  TJ-l 

“Pi-i - <*1 

The  diagonal  block  Rp,_i.Pi_i  may  be  written  as 

R-Pi-i.p,-i  =  T  +  Ap,_i,p,_i 


(3.33) 


(3,34) 


where  is  the  difference  between  Rpj_i^,_i  and  the  Toeplitz  approximation 

T. 

Substituting  (3.34)  into  (3.33)  modifies  the  recursion  to  give 

“l‘’  =  Gr'h, -h, (335) 

and 

-  T-'hfai"-') .  (3.36) 

The  final  expressions  resulting  from  the  successive  pau’titioning  are  then 

“Pi-i  =  ^p^-iK  p,-i  —  hpi-i  ,  (3.37) 

and 

=  -Rx.,hp,.,a<p*:V  •  (3.38) 

where  ap,_i  =  ao,Gp,_i  =  Ro,o,  and  hpj_i  =  Ro.i.  The  diagonals  Rjj  may  be 
rewritten  as 

Rjj  =  T  +  Ajj  ,  (3.39) 
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Now,  substituting  (3.39)  into  (3.38)  results  in 

-  hn-,  .  (3.40) 

and 

^  .  (3.41) 

Finally,  combining  the  succesive  iterative  equations  results  in  our  solution  to  the 
normal  equations  using  the  forward  iterative  algorithm  for  the  QP  case,  which  takes 
the  form 

a?‘'  =  .?">  -  T- AoX"*"'  -  T-  s'  (3.42) 

t=l 

af*l  =  -T-' AyjaJ<‘-*>  -  T->  s'  (3.43) 

imO 

i*3 

where  =  T“^Rj,oaJ[^®^  and  j  =  1,2, . . .  ,Pi  —  1.  The  forward 

iteration  will  provide  parameters  for  the  first  quadrant  spectral  estimate. 

2.  Backward  Iteration 

A  filter  mask  may  have  quarter-plane  support  in  any  two  quadrants.  A  sec¬ 
ond  quadrant  estimate  exists  for  the  quadrant  adjacent  to  either  side  of  the  quadrant 
designated  to  be  the  first  quadrant.  The  Hermitian  symmetry  and  Toeplitz-block- 
Toeplitz  nature  of  the  autocorrelation  matrix  causes  the  estimates  produced  from 
diagonally  adjacent  quadrants  to  be  identical  [Ref.  9]. 

The  second  quadrant  AR  parameters  6(t,  j)  may  be  estimated  from  the  nor¬ 
mal  equations  using  the  backward  iterative  Toeplitz  approximation.  When  6(0, 0)  =  1 
the  second  quadrant  normal  equation  is  given  by 


Ro.o  Ro.i  Ro,2  •"  Ro^i-i 

'bp,_i  ■ 

0 

R-1,0  R-i.i  R'l.a  ••• 

bpi-2 

0 

R-2,0  B<2,1  R-2,2  •  •  R2,P,_1 

bpi-3 

= 

0 

•  •  •  •  • 

.R-Pi-1,0  Rpi-l.l  R-Pi-1,2  •••  B-P,-l.Pi-l. 

.  bo  . 
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where  =  [100,...,  0]^  and 

bk  =  (!/<’’) Wt,  0),  6(«:,  1),  6(t,  2), . . . ,  Ki,  Pi  -  1)1 . 


It  may  be  noted  that  the  second  quadrant  autocorrelation  matrix  is  identical  to  that 
of  the  first  quadrant.  However,  the  parameter  estimates  b{i,j)  aie  not  generally  the 
same.  The  estimation  of  the  backward  parameters  is  computed  iteratively  in  the  same 
manner  as  the  forward  parameters,  but,  the  backward  method  differs  in  the  way  that 
the  normal  equation  is  partitioned.  In  this  case  the  partitioning  begins  at  the  upper 
left  diagonal  block  and  continues  until  the  G  matrix  consists  of  only  the  lower  right 
diagonal  block.  The  Toeplitz  approximation  T  is  identical  to  that  of  first  quadrant 
support. 

The  combined  backward  iterative  equations,  which  £u^e  similar  to  those  for 
the  forward  iteration,  can  be  written  as 


b'*'  =  bf  -  T- An..,P.-,b<‘->  -  T-  t  Rp..,^...ibn.i 


(3.45) 


bf  =  -  T->  £  (3-46) 


where  bo**  =  bj®*  =  T"*Rr,_ij'bo®^  and  j  =  1,2, .. .  ,Pi  —  1.  The  first  and 

second  quadrant  parameters  are  used  to  find  the  combined-quadrant  spectral  estimate, 
which  is  discussed  in  Chapter  V. 


3.  Iterative  Method  for  NSHP  Support 

The  iterative  Toeplitz  approximation  method  may  be  used  to  estimate  the 
parameters  of  arrays  with  nonsymmetric  half-plane  support.  As  shown  previously  the 
NSHP  autocorrelation  matrix  has  an  asymmetric  structure  which  can  be  expressed 
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(3.47) 


Ri,o 

R-1,1 

R-1.2 

R-i,Pi-i 

R2,0 

R-2,1 

R’2,2 

R-2,Pi-1 

.Rpi_i,o 

R-Pi-1.1 

R.Pi-1,2 

•  * 

Rpj_i,Pj_ 

The  dimensions  of  the  top  most  and  left  most  blocks  axe  reduced  by  rows  and 
columns,  respectively.  As  in  the  previous  caaes  the  iterative  process  begins  with 
forming  the  Toeplitz  approximation  of  Ravy  In  the  NSHP  case,  however,  the  upper 
left  main  diagonal  block  is  not  included  in  the  average  becuaae  it  does  not  have 
the  same  dimension  as  the  other  main  diagonal  blocks.  Thus,  it  is  necessary  to  use 
the  elements  of  the  upper  left  msdn  diagonal  block  to  create  a  separate  Toeplitz 
approximation  designated  as  T. 

The  normsJ  equation  is  successively  partitioned  as  described  in  the  previous 
sections.  At  each  stage  of  the  partitioning,  the  blocks  along  the  main  diagonal  are 
expressed  as  the  sum  of  the  Toeplitz  approximation  and  a  difference  matrix.  The 
iterative  sununations  are  given  by 


- 1-'  s’ 

«=1 

(3.48) 

Pj-1 

_  T-1 

fsO 

(3.49) 

where  and  j  =  1, 2, . . . ,  A  -  1. 


C.  SUMMARY 

It  is  clear  that  iterative  methods  have  a  computational  advantage  over  the  di¬ 
rect  inversion  method.  This  has  been  shown  for  the  Toeplitz-block-Toeplitz  case  and 
extended  to  the  case  where  Toeplitz  approximation  is  used  to  compensate  for  the 
non-Toeplitz  nature  of  autocorrelation  matrices  formed  using  the  covariance  method. 


Generally,  the  matrix  that  must  be  inverted  using  the  iterative  method  has  the  di¬ 
mension  of  a  single  block  in  the  block  matrix.  This  corresponds  to  the  filter  order  or 
mask  size.  The  mask  size  will  vary  with  the  order  necessaxy  for  accurate  AR  mod¬ 
eling  of  the  sampled  data,  but  filter  masks  that  range  from  3  x  3  to  6  x  6  will  be 
sufficient  for  most  applications.  Although  other  methods  may  exhibit  faster  conver¬ 
gence  to  the  true  parameters,  they  do  so  with  greater  computational  complexity.  The 
iterative  methods  provide  an  approximation  much  sooner.  Additionally,  the  iterative 
methods  are  guaranteed  to  converge  when  the  block  matrix  is  symmetric  and  positive 
definite  [Ref.  2].  This  property  is  importamt  when  extending  the  iterative  Toeplitz 
approximation  algorithm  to  the  data-adaptive  case  which  is  introduced  in  the  next 
chapter. 
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IV.  DATA- ADAPTIVE  ITERATIVE  TOEPLITZ 
APPROXIMATION  ALGORITHM 


A.  OVERVIEW 

The  necessity  for  adaptive  filtering  techniques  arises  when  it  is  desired  to  process 
signals  that  result  from  an  environment  of  unknown  statistics  [Ref.  10].  There  exists  a 
broad  class  of  problems  that  fall  into  this  category,  these  include  such  diverse  fields  as 
sonar,  radar,  image  processing,  seismology  and  communications.  In  general,  adaptive 
filters  provide  a  significant  improvement  in  performance  over  fixed  coefficient  filters 
designed  to  operate  on  data  with  unknown  statistics.  Additionally,  the  use  of  adaptive 
filters  makes  possible  new  signal  processing  capabilties  that  would  not  be  available 
otherwise. 

One  distinct  advantage  associated  with  adaptive  filters  pertains  to  their  ability 
process  data  on*line.  This  is  desirable  for  many  applications  such  as  autoregressive 
spectrum  analysis,  detection  of  a  signal  in  noise,  adaptive  noise  cancellation  and  Une 
enhancement. 

Two-dimensional  adaptive  filters  are  used  to  process  data  obtained  from  spa¬ 
tially  variant  data  arrays.  The  most  successful  algorithms  currently  implemented  for 
this  purpose  include  the  2-D  least  mean  square  (LMS)  and  recursive  least  squares 
(RLS)  algorithms.  Both  of  these  algorithms  have  been  derived  from  the  2-D  Wiener 
filter  and  method  of  least  squares  [Ref.  11].  With  this  in  mind,  it  follows  that  the 
iterative  method  for  fixed  data  may  be  successfully  extended  to  operate  on  spatially 
variant  arrays.  The  remainder  of  this  chapter  describes  the  development  of  the  data- 
adaptive  iterative  Toeplitz  approximation  algorithms  and  some  considerations  that 
apply  when  using  this  algorithm. 
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B.  THE  DATA- ADAPTIVE  ALGORITHM 


The  data-adaptive  Toeplitz  approximation  algorithm  is  based  on  the  same  suc- 
cesive  partitioning  described  for  the  fixed  data  case.  The  algorithm  becomes  adaptive 
in  nature  when  the  autocorrelation  matrix  R„  is  updated  for  each  shift  n  of  the  filter 
mask.  This  yields  a  continuously  updated  set  of  AR  parameter  estimates.  Compu¬ 
tation  of  AR  parameter  estimates  begins  with  the  first  iteration  [Ref.  12].  As  in  the 
fixed  data  case,  only  the  approximated  Toeplitz  matrix  T  is  inverted.  Computation 
time  is  shortened  further  by  the  fact  that  R„  is  computed  only  from  data  covered 
by  the  filter  mask,  which  amounts  to  taking  the  outer  product  of  a  PiPj  x  1  vector 
where  PiPj  is  equal  to  the  number  of  filter  coefficients.  The  data-adaptive  form  of 
the  combined  iterative  equations  can  be  written  as 
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where  aj™  =  a„"*"'  =  T„-*R„,,oa?°’ “d  J  =  1,2,. . .  ,F, -1.  These  can  be 

compajed  to  (3.42)  and  (3.43)  for  the  fixed  data  case.  The  initial  parameter  estimates 
sto*'  Jire  determined  for  the  initial  point  being  estimated  and  all  subsequent  parameters 
a"  aie  a  function  of  the  parameters  computed  at  the  previous  mask  positions.  The 
backward  equations,  used  for  the  second  quadrant  estimates,  cem  be  written  as 
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where  =  T;^fW,  bof  =  To-'RoP.-ijbi®)  and  ;  =  1,2, . . . , Px  -  1. 

It  should  be  noted  that  the  only  difference  between  the  data-adaptive  equations 
and  the  fixed  data  equations  is  the  the  index  of  iteration  which  is  n  instead  of  k 
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and  the  manner  in  which  the  autocorrelation  matrix  is  formulated.  In  the  data 
adaptive  case,  the  iteration  proceeds  with  each  shift  of  the  filter  mask  or  update  of 
the  autocorrelation  matrix  which  has  the  recursive  form 

R„  =  R„_1  +  -Knxl .  (4.5) 

The  P1P2  X  1  vector  Xn  is  obtained  by  arranging  the  2-D  data  covered  by  the  filter 
mask  at  each  shift  n  =  (ni,  nj)  in  a  vector  form.  The  diagonal  blocks  of  Rn  are  used 
to  form  Tn  in  the  same  manner  as  the  fixed  data  case. 

As  with  the  RLS  algorithm  and  other  recursive  implementations  of  the  method 
of  least  squares,  it  is  necessary  to  introduce  a  weighting  or  forgetting  factor  when 
computing  R„  [Ref.  10].  The  forgetting  factor  is  a  scalar  value  that  may  be  designated 
as  ^{n)  where  n  is  the  iteration  number  of  the  point  being  estimated.  The  weighting 
factor  ^{n)  has  the  property  that 

0  < /?(n)  <  1,  n  =  l,2,...,p,  (4.6) 

where  p  is  the  total  number  of  points  being  estimated.  The  purpose  of  ^{n)  is  to 
ensure  that  data  in  the  distant  past  is  ^forgotten’  which  will  make  it  possible  for 
the  filter  to  adapt  to  statistical  variations  of  the  observed  data  when  operating  in  a 
non- homogeneous  environment.  A  commonly  used  form  of  the  weighting  factor  is  the 
exponential  weighting  factor  which  is  defined  as 

/3(n)  =  A"'”,  n  =  l,2,...,p  ,  (4.7) 

where  A  is  positive  constant  with  a  value  less  than  1.  The  quantity  1/(1  —  A)  can  be 
considered  as  2m  approximate  measure  of  the  memory  of  the  algorithm.  When  A  =  1 
we  have  the  ordinary  method  of  least  squares  which  corresponds  to  infinite  memory. 
By  this  reasoning,  we  may  introduce  a  method  of  exponentially  weighted  least  squares, 
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where  y9(n)  is  used  as  a  weighting  factor  in  the  expression  of  the  performance  measure 
[Ref.  10] 

«")  =  E  A'-"|e(n)l’  ,  (4.8) 

n=l 

where  e(n)  is  the  error  defined  by  (2.12)  at  mask  position  (ni,n2). 

To  see  how  the  forgetting  factor  is  introduced  in  the  recursion  we  must  refer  to 
the  original  definition  of  the  normal  equation 

Ra  =  £  ,  (4.9) 

where  a  is  the  vector  of  optimum  parameters  which  results  in  the  minimum  value  of 
the  performance  measure  ^(n).  The  PiPj  x  P1P2  autocorrelation  matrix  in  this  case 
is  defined  by 

Rn  =  •  (4.10) 

n=l 

We  may  now  isolate  the  term  of  (4.10)  that  corresponds  to  n  =  p  from  the  rest  of  the 

summation  on  the  right  side  of  the  expression  to  obtain 

V-i 

Rn  =  A  51  +x„x^.  (4.11) 

n=sl 

The  term  inside  the  brackets  on  the  right  side  of  (4.11)  is  actually  the  sum  of  the  pre¬ 
viously  computed  autocorrelation  matrices  R(n  —  1),  therefore  we  have  the  recursive 
relationship 

Rt,  =  AR„_i  -f  x„x^  .  (4.12) 

The  affect  of  ^{n)  is  apparent  when  (4.12)  is  expanded  to  give 

Ri  =  xixf 

R2  =  Axixf  -I-  X2X2 

R3  =  A^xixf  +  AX2X2  +  xsxj 

Rn  =  A"~^Xixf -1- A""^X2X2  +  . . .  +  X„x^  (4.13) 
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which  clearly  shows  the  decreasing  contribution  made  by  past  data  vectors  when 
A  <  1. 

The  choice  of  a  value  for  A  is  dependent  on  the  statistical  nature  of  the  data 
being  processed.  Generally,  if  the  data  is  random  or  uncorrelated  it  is  undesirable  to 
use  the  past  data  to  compute  the  current  estimate,  therefore  a  small  value  for  A  would 
be  the  appropriate  choice.  Conversely,  a  value  of  A  close  to  1  would  be  desirable  for 
highly  correlated  data. 

Since  the  statistical  nature  of  the  data  is  often  not  known  a  priori  when  using 
data-adaptive  hlters,  it  is  difficult  to  choose  an  optimum  value  for  the  forgetting 
factor.  A  further  modification  of  the  recursive  equation  (4.12)  has  been  implemented 
for  the  data-adaptive  iterative  algorithm  to  address  this  problem.  This  modification 
involves  the  application  of  a  weighting  factor  to  the  current  update  of  Rn.  This 
weighting  factor  depends  on  the  index  of  iteration  n  such  that  its  value  approaches 
1  as  n  — ►  oo.  This  biasing  scheme  causes  the  influence  of  the  present  data  on  the 
parameter  estimates  to  increase  with  the  number  of  iterations.  The  modified  recursion 
is  given  by 

R„  =  AR„_i  +  XnX^  .  (4.14) 

The  weighting  factor  on  the  second  term  in  (4.14)  is  added  to  provide  a  more  gradual 
shift  of  influence  to  recent  updates  than  wotild  be  possible  with  only  the  A  term 
present.  This  factor  is  not  as  important  when  a  highly  correlated  signal  is  processed. 

C.  SUMMARY 

The  adaptive  Toeplitz  approximation  algorithm  can  be  used  to  estimate  pauram- 
eters  for  the  same  regions  of  support  as  in  the  fixed  data  case.  Additionally,  various 
schemes  for  pre-windowing  or  post-windowing  the  data  may  be  employed  to  improve 
the  estimates.  Other  factors  affecting  the  performance  of  the  algorithm  include  the 
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directional  scheme  for  moving  the  filter  over  the  data  array  and  the  criteria  used  to 
determine  convergence  of  the  parameter  estimates.  This  algorithm  has  been  used  for 
spectral  estimation,  image  noise  cancellation,  and  line  enhancement.  Experimental 
results  involving  these  applications  are  presented  in  the  next  two  chapters. 
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V.  2-D  AR  SPECTRAL  ESTIMATION 


A.  OVERVIEW 

One  application  where  the  use  of  two-dimensional  autoregressive  modeling  has 
proven  to  be  particularly  advantageous  is  that  of  spectr2J  estimation.  The  2-D  AR 
model  was  derived  in  Chapter  II.  This  model  is  used  to  find  the  2-D  AR  power 
spectral  estimate  which  is  given  by  [Ref.  4] 

=  \H{uh,u)2)\^Ptv  ,  (5.1) 


where  P^,  is  the  power  spectral  density  of  the  input  and  is  the  transfer 

function  of  the  2-D  AR  model.  The  input  is  white  noise  with  a  constant  power 
spectrum  of  amplitude  al,.  Therefore,  we  can  write  (5.1)  as 


(5.2) 


The  key  to  finding  the  2-D  AR  power  spectral  estimate  is  determining  the  the  pa¬ 
rameters  a{ki,k2).  The  next  section  provides  results  of  experimentation  with  the 
data-adaptive  iterative  algorithm  as  used  to  find  the  parameter  estimates  of  a  2-D 
signal  consisting  of  a  sum  of  sinusoids  in  additive  white  noise.  Results  using  the  di¬ 
rect  inversion  method  and  fixed  data  iterative  Toeplitz  approximation  are  provided 
for  comparison. 


B.  EXPERIMENTAL  RESULTS 

The  performance  of  the  algorithms  discussed  in  the  previous  chapters  is  com¬ 
pared  by  examining  the  estimated  spectral  densities  computed  from  the  parameters 
produced  by  each  method.  Each  algorithm  has  been  used  to  estimate  the  parameters 
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of  the  sinusoidal  input  array  generated  by 


z(ni,  n2)  =  Ai  cos(27r/nni  +  27r/i2n2) 

+>l2  cos(27r/2ini  +  2irf22n2)  +  w{ni,n2)  (5.3) 

where  amplitudes  Ai  =  A2  =  y/2  and  /,j  axe  normalized  frequencies  in  the  range 
(0  ^  fij  ^  0-5).  The  frequencies  chosen  are  fn  =  /12  =  0.167  and  /21  =  f22  =  0.333. 
Additioncilly,  the  data-adaptive  algorithm  is  tested  with  off-set  frequencies  to  evaluate 
its  peformzmce  for  that  case.  The  noise  term  tw(ni,n2)  is  zero  me2Ln  and  gaussian  with 
the  variance  scaled  give  a  desired  signal-to-noise  ratio  (SNR).  The  SNR  (in  dB) 
is  defined  as  [Ref.  4] 

SNR=101og,„£;^,  (5.4) 

where  N  is  the  number  of  sinusoids  present  and  A  is  the  amplitude  term.  The  SNR 
is  varied  by  holding  A,-  constant  and  adjusting  the  value  of  o-^,.  The  algorithm’s 
performance  using  a  3  x  3  filter  mask  is  evaluated  for  various  sizes  of  input  arrays 
and  the  common  regions  of  support,  including  combined-quadrant  (CQ).  The  results 
obtained  using  the  data-adaptive  algorithm  are  provided  for  the  pre-windowed  eind 
covariance  methods.  The  data-adaptive  algorithm  is  tested  with  SNR’s  of  10  dB  and 
0  dB.  It  should  be  noted  that  the  surface  plots  in  all  figures  have  been  rotated  90 
degrees  to  show  the  separation  of  the  spectral  peaks  more  clearly.  Contour  plots  axe 
provided  to  better  judge  the  accuracy  of  the  estimates.  The  axes  of  the  contour  plots 
range  from  0-30  which  is  taken  from  a  32  x  32  frequency  domain  array  representing  a 
quadrant  of  a  64  x  64  point  2-D  FFT  output.  The  32  x  32  array  covers  the  frequency 
range  from  0  to  tt. 

1.  Estimates  by  Direct  Inversion 

The  first  case  examined  is  that  involving  direct  inversion  of  the  autocorre¬ 
lation  matrix  to  find  the  least  squares  AR  parameters  from  the  normal  equations.  An 
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8x8  input  array  was  used  with  1*‘  quadrant  quarter-plane  support.  The  resulting 
spectral  density  estimate  is  shown  in  Figure  5.1.  This  plot  shows  spectral  peaks  at 


Figure  5.1:  Spectral  estimate  using  direct  matrix  inversion;  SNRslO  dB. 

21.3  and  10.7  on  the  Fl  and  F2  axis  respectively,  which  correspond  to  the  normalized 
frquencies  /u  =  /u  =  0.167  and  /ji  =  /m  =  0.333  of  the  input  data.  The  true 
locations  of  the  frequencies  are  indicated  by  crosses. 

2.  Spectral  Estimates  using  the  Iterative  Method  for  Fixed  Data 
a.  Quarter>Plane  Support 

Figure  5.2  shows  the  spectral  density  of  3;(ni,n2)  as  estimated  using 
the  fixed  data  iterative  method  with  first  quadrant  QP  support.  As  in  the  previous 
case  the  frequencies  estimated  are  very  close  to  the  true  frequencies.  It  is  interesting 
to  note  however,  that  contrary  to  what  might  be  expected,  the  quality  of  the  estimate 
appears  to  degrade  as  the  number  of  iterations  is  increased.  The  estimate  produced 
after  one  iteration  is  clearly  superior  to  the  estimate  obtained  after  ten  iterations. 
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This  phenomenon  is  consistent  with  that  experienced  in  the  1-D  case  where  the 
iterative  solution  would  begin  to  diverge  after  a  certain  number  of  iterations  in  cases 
with  small  data  sets  [Ref.  12]. 


I  liMtklMI 
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Figure  5.2:  Fixed  Data,  1st  Quadrant  QP  Support,  SNR=10  dB. 

Results  using  second  quadrant  QP  support  proved  to  be  generally 
better  for  this  placement  of  the  sinusoids  than  those  obtained  from  the  First  quadrant. 
An  example  is  provided  in  Figure  5.3. 

b.  Combined  Quadrant  Support 

Spectral  density  estimates  using  QP  support  demonstrate  a  tendency 
to  distort  elliptically.  This  distortion  is  evident  in  the  first  quadrant  QP  case  shown 
in  Figure  5.2.  A  method  to  compensate  for  this  elongation  of  the  spectral  peaks 


31 


involves  combining  the  first  and  second  quadrant  estimates  [Ref.  13]  to  give 


The  plots  in  Figure  5.4  indicate  that  this  method  yields  accurate  estimates.  However, 
despite  some  improvement,  the  elliptical  elongation  at  the  base  of  the  spectral  peaks 
resulting  from  the  influence  of  the  first  quadrant  estimate  is  still  quite  noticeable. 


Figure  5.4:  Fixed  Data,  Combined  Quadrant  Support,  SNR=10  dB. 

3.  Spectral  Estimates  using  the  Data-Adaptive  Iterative  Method 
a.  Pre- Windowed  Input  Arrays 

The  input  arrays  were  pre-windowed  by  adding  Pi  —  1  rows  and 
rol limns  of  zeros  to  the  top  and  left  side  of  the  array,  where  Pi  is  the  dimension 
of  the  filter  mask.  The  purpose  of  pre-windowing  the  data  is  to  compute  the  pa¬ 
rameter  estimates  using  all  the  input  data.  Figure  5.5  shows  the  spectral  estimates 
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for  1*‘,  and  CQ  support  with  an  8  x  8  input  array.  In  Figure  5.6  we  present  the 


Figure  5.5:  8  x  8  input  array  (pre- windowed),  SNR=10  dB. 


spectral  estimates  when  only  16  data  points  (4x4  array)  were  available.  The  spectral 
estimates  obtained  from  the  8x8  input  array  were  quite  acceptable,  but  some  bias 
is  evident  in  the  first  quadrant  estimate  which  is  carried  over  to  the  CQ  estimate. 
The  estimates  obtuned  from  the  4x4  array  demonstrate  the  value  of  computing  the 
CQ  estimate.  In  this  case  the  true  signal  frequencies  were  unrecognizable  until  the 
combined  quadrant  estimate  was  obtained. 

A  final  test  of  the  pre-windowed  algorithm  consisted  of  running  the 
filter  mask  only  half  way  through  an  8  x  8  input.  The  estimates  are  computed  from 
the  upper  4x8  section  of  the  observed  data.  The  result  of  this  experiment  is  given  in 
Figure  5.7.  These  results  are  similar  to  that  of  the  4x4  case  and  serve  to  reinforce  the 
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Figure  5.7:  4x8  portion  of  the  input  array  (pre- windowed),  SNR=:10  dB. 


importance  of  the  Combined  Quadrant  technique  for  2-D  AR  spectrum  estimation, 
b.  Estimates  using  the  Covariance  Method 

As  explained  in  section  B  of  Chapter  II  the  covariance  method  uses 
only  the  observed  data.  As  a  result,  Pi  —  1  rows  and  columns  of  the  input  array  are 
not  available  when  computing  the  parameter  estimates.  However,  the  experimental 
results  obtained  in  this  work  indicate  that  the  covariance  method  achieves  better 
spectral  peak  resolution.  Additionally,  less  bias  is  observed  in  the  spectral  estimates. 
This  characteristic  of  the  covariance  method  is  evident  in  Fig  5.8.  In  this  case  the 


Figure  6.8:  8x8  input  array  (Covariance  Method),  SNR=10  dB. 

best  spectral  estimate  is  clearly  provided  by  the  Combined  Quadrant  method.  The 
bias  present  in  the  pre-windowed  first  quadrant  support  example  has  been  completely 
eliminated  by  using  the  covariance  method. 
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The  next  example  demonstrates  a  unique  characteristic  of  the  data- 
adaptive  iterative  Toeplitz  approximation  algorithm.  In  this  example  the  spectral 
estimate  is  computed  using  a  3  x  3  filter  mask  with  a  4  x  4  input  array.  This  means 
that  only  four  elements  of  the  input  data  are  available  to  the  filter  to  determine 
the  parameter  estimates.  This  result  is  of  particular  interest  because  the  correlation 
matrix  is  actually  rank  deficient,  and  the  problem  cannot  be  solved  by  direct  inversion. 
Nevertheless,  the  iterative  Toeplitz  approximation  method  does  provide  a  means  to 
compute  the  estimate.  The  results  are  given  in  Figure  5.9.  Although  the  spectral 
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Figure  5.9:  4x4  input  array  (Covariance  Method),  SNR=10  dB. 

estimates  for  the  first  quadrant  and  Combined  Quadrants  are  somewhat  weak,  the 
second  quadrant  estimate  exhibits  remarkable  resolution  and  accuracy. 

As  a  further  test  of  the  algorithm,  a  signal  with  offset  frequencies 
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fn  =  0.250,  /i2  =  0.125,  /21  =  0.300,  and  /22  =  0.400  is  generated  and  processed. 
The  spectral  estimates  resulting  from  the  computed  model  parameters  are  plotted 
in  Figure  5.10.  For  this  caise,  the  best  spectral  estimate  was  obtained  using  first 
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Figure  5.10:  8x8  input  array,  offset  frequencies,  (Covariance  Method), 
SNR=10  dB. 

quadrant  suppport. 

The  spectral  estimates  for  a  signal-to-noise  ratio  of  0  dB  is  ploted 
in  Figure  5.11.  The  accuracy  of  the  estimates  is  noticeably  degraded  in  this  case. 
However,  a  reasonable  indication  of  the  signaj  frequency  content  can  be  seen  when 
the  combined-quadrant  method  is  used. 

The  final  ceise  using  NSHP  support  is  included  for  completeness.  The 
filter  mask  had  dimensions  Pi  =  P2  =  i  and  L2  =  —2.  The  parameters  were  estimated 
for  a  16  X  16  input  array.  The  NSHP  spectral  estimate  is  plotted  in  Figure  5.12.  The 
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resulting  spectral  estimate  shown  in  Figure  5.12  is  quite  accurate,  but  this  comes 
at  the  cost  of  greater  complexity  for  the  filter.  NSHP  support  may  be  used  as  an 
alternative  to  combined-quadrant  support  in  some  cases. 

C.  SUMMARY 

The  data-adaptive  iterative  Toeplitz  approximation  cdgorithm  has  proven  to 
be  a  viable  means  for  2-D  AR  spectrum  estimation.  Numerous  factors  must  be 
considered  when  using  this  algorithm.  In  particular,  the  filter  mask  size  should  be 
chosen  to  correspond  to  the  model  order  of  the  signal  being  processed.  Although 
our  overall  experience  has  found  that  that  CQ  support  produced  the  best  results, 
no  region  of  support  seemed  to  have  a  cleau-  advantage  over  the  others  for  all  cases. 
This  indicates  that  the  optimum  region  of  support  must  be  evtiluated  for  the  specific 
application  that  the  aJgorithm  is  implemented  for. 
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VI.  RESTORATION  OF  IMAGES 

A.  OVERVIEW 

The  characteristics  of  typical  images  vary  greatly  from  one  region  of  the  image 
to  the  next.  For  example,  a  region  of  an  image  contzuning  a  crowd  or  a  building  may 
have  detailed  variations  in  intensity,  while  another  region  representing  the  sky  in  the 
backgroimd  will  be  essentially  uniform  in  intensity.  Two-dimensional  data-adaptive 
filtering  is  particularly  suited  for  processing  data  of  this  nature.  The  idea  of  adapting 
the  processing  to  the  local  characteristics  of  an  image  is  sidvantageous  for  many  image 
processing  applications,  including  image  enhancement  and  restoration. 

There  are  two  approaches  to  adaptive  signal  processing.  The  first  approach 
involves  adapting  the  Alter  output  for  each  pixel.  This  is  known  as  pixel-by-pixel 
processing  [Ref.  4].  In  this  scheme,  the  processing  method  is  based  on  the  local 
characteristics  of  the  image,  degradation,  and  any  other  pertinent  information  con¬ 
tained  in  the  pixel’s  neighborhood  region.  This  approach  offers  the  greatest  degree 
of  flexiblity  to  the  adaptive  process,  but  has  the  highest  computational  complexity. 

When  a  more  computationally  efficient  method  is  desired,  a  second  approeich 
known  as  suhimage-hy-snbimage  or  block-by-block  processing  can  be  used  [Ref.  4].  In 
this  approach,  the  image  array  is  divided  into  subimages  and  space-invariant  tech¬ 
niques  are  used  to  process  the  data.  This  method  can  result  in  some  discontinuities 
being  present  in  the  processed  image.  The  si25e  of  the  subimages  directly  affects  the 
quality  of  the  output  image.  The  block-by-block  method  provides  a  faster  means  of  ^ 

processing,  however,  and  may  be  desirable  in  some  applications. 

In  the  next  two  sections  of  this  chapter  the  data-adaptive  Toeplitz  approxima- 
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tion  algorithm  is  implemented  in  image  noise  cancellation  and  line  enhancer  modes. 
For  both  of  these  applications  the  algorithm  has  been  modified  to  perform  pixel-by¬ 
pixel  processing. 

B.  IMAGE  NOISE  CANCELLATION 

The  usual  method  of  finding  the  estimate  of  a  signal  in  noise  is  to  pass  the 
corrupted  signal  through  a  system  that  serves  to  suppress  the  noise  while  leaving 
the  desired  signal  relatively  unchanged.  The  noise  canceler  developed  by  Widrow 
(Ref.  14]  is  an  example  of  such  a  system.  Adaptive  noise  cancellation  is  a  variation 
of  optimal  filtering  that  can  be  used  for  image  restoration  to  some  advantage.  In 
particular,  when  a  recieved  image  is  obscured  by  noise  or  another  image  transmitted 
from  another  location  as  shown  in  Figure  6.1,  an  adaptive  noise  canceler  can  be 
employed  to  extract  the  desired  image  from  the  composite  signal.  The  noise  canceler 


unwanted  transmission 


Figure  6.1:  A  possible  scenario  for  application  of  a  noise  canceler  as  used 
for  image  restoration. 

makes  use  of  an  auxiliary  or  reference  input  derived  from  one  or  more  sensors  located 
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at  points  in  the  noise  field  where  the  desired  signal  is  weak  or  undetectable.  This 
input  is  filtered  and  subtracted  from  the  primary  input  containing  both  desired  and 
unwanted  signals.  A  block  diagram  of  the  noise  canceler  is  shown  in  Figure  6.2.  The 
reference  signal  and  the  undesired  part  of  the  primary  input  are  correlated.  For  this 
reason,  the  unwanted  signal  or  noise  is  eliminated  or  attenuated  by  cancellation. 


Figure  6.2:  Adaptive  Noise  Canceler  Block  Diagram. 


To  evaluate  the  performance  of  the  data-adaptive  iterative  Toeplitz  approxi¬ 
mation  algorithm  used  in  a  noise  canceler  configuration,  a  256  x  256  composite  im¬ 
age  array  was  created  from  two  distinct  images.  This  corrupted  image  is  shown  in 
Figure  6.3.  The  algorithm  described  in  Chapter  IV  was  modified  to  include  a  cross¬ 
correlation  term  given  by  ^ 

r„  =  Ar„_i  +  (^-1)  > 

where  <f„  is  the  element  of  the  composite  signal  corresponding  to  the  element  of 
the  reference  signal  being  processed  by  the  filter.  This  cross-correlation  term  is  used 
to  form  the  initial  parameter  estimate  for  each  pixel  being  processed. 
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Figure  6.3:  Corrupted  image  array  provided  as  primary  input  to  noise 
canceler. 

The  resulting  error  e(ni,n2)  computed  at  the  output  of  the  system  becomes  the 
restored  image.  This  relationship  is  given  by 

e(ni,n2)  =  [d(ni,n2)  +  u(ni,n2)]  -u(ni,n2)  ,  (6.2) 

where  d  is  the  desired  image,  it  is  the  unwanted  image  or  noise  and  u  is  the  filtered 
estimate  of  u. 

In  the  simulation  the  noise  signal  added  to  the  desired  image  is  scaled  and 
lowpass  filtered  to  simulate  the  effect  of  having  the  primeiry  and  reference  sensors 
physically  separated.  The  output  of  the  noise  canceler  is  shown  in  Figure  6.4.  The 
unwanted  image  h2is  been  noticeably  attenuated  and  the  desired  image  is  clearly 
visible.  The  original  images  are  provided  in  Figure  6.5  as  a  reference  for  evaluation 
of  the  algorithm’s  perform2uice. 

C.  ADAPTIVE  LINE  ENHANCER 

A  special  case  of  adaptive  noise  canceling  occurs  when  only  the  signal  that  is 
corrupted  by  noise  is  available.  In  this  case  the  recieved  signal  is  delayed  and  used  as 
the  reference  signal.  The  reference  signal  may  be  expressed  as  y(ni,  n2)  =  a:(ni  —  6, 712), 
where  6  is  the  amount  of  spatial  delay  used.  The  block  diagram  for  this  system  is 
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shown  in  Figure  6.6.  The  main  function  of  the  delay  parameter  S  is  to  remove  any 


Figure  6.6:  The  adaptive  line  enhancer  block  diagram. 

correlation  that  may  exist  between  the  noise  component  in  the  original  input  signal 
and  the  noise  component  of  the  delayed  adaptive  filter  input  [Ref,  10].  The  filter  will 
respond  by  cancelling  any  components  of  the  main  signal  x(ni,n2)  that  are  in  any 
way  correlated  with  the  secondary  signal  y(ni,n2)  [Ref.  11].  The  remaining  noise 
component  at  the  output  of  the  filter  is  then  subtracted  from  main  signal  resulting 
in  the  removal  of  the  additive  noise  in  the  signal.  In  general,  the  delay  6  should 
be  chosen  such  that  it  is  approximately  half  the  length  of  the  correlation  sequence 
[Ref.  11].  The  input  image  and  output  of  the  adaptive  line  enhancer  are  shown  in 
Figure  6.7. 

D.  SUMMARY 

In  this  chapter,  it  has  been  demonstrated  that  the  data  adaptive  Toeplitz  ap¬ 
proximation  algorithm  may  be  used  in  image  processing  applications.  While  pixel-by¬ 
pixel  processing  was  used  in  the  above  examples,  it  may  be  desirable  to  combine  this 
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Figure  6.7:  Input  and  output  of  the  data-adaptive  line  enhancer. 


approach  with  block-by-block  processing  to  obtsun  better  results.  This  would  entail 
dividing  the  image  into  smaller  blocks  before  processing  and  then  recombining  the 
processed  blocks.  Additionally,  reprocessing  of  the  sirray  estimates  may  be  a  viable 
means  of  improving  image  quality. 
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VII.  CONCLUSIONS 


The  data-adaptive  iterative  Toeplitz  approximation  algorithm  is  readily  im¬ 
plemented  from  the  iterative  methods  previously  used  for  fixed  data-arrays.  The 
experiment2il  results  indicate  that  this  algorithm  is  well-suited  for  a  wide  variety  of 
2-D  signal  processing  applications.  Practically  speaking,  the  iterative  method  can 
be  successfully  implemented  for  any  application  where  the  use  of  Wiener  filter  bcised 
adaptive  filters  has  been  successful.  While  more  extensive  evaluation  is  necessary 
before  this  method  can  be  said  to  have  a  clear  advantage  over  other  methods,  there 
is  enough  evidence  of  its  capabilities  to  warrant  further  investigation  in  this  area. 

A.  PERFORMANCE  EVALUATION  SUMMARY 

As  stated  earlier,  autoregressive  model  parameters  may  used  to  estimate  the 
spectral  content  of  2-D  arrays  with  high  resolution.  In  experiments  the  performeince 
of  the  data-adaptive  iterative  algorithm  matched  that  of  the  direct  inversion  method, 
2Lnd  in  some  cases,  improved  upon  the  performance  of  the  fixed  data  iterative  method. 
This  improvement  may  be  due  to  the  way  in  which  the  iteration  is  carried  out  in  the 
adaptive  caae  versus  the  fixed  case.  In  the  fixed  case  the  iteration  is  carried  out  using 
a  correlation  matrix  formed  from  all  the  data  at  once.  In  the  data-adaptive  case  the 
correlation  matrix  is  formed  recursively  with  new  data  incorporated  at  each  itera¬ 
tion.  Of  particular  note,  W2is  the  ability  of  the  data-adaptive  eilgorithm  to  estimate 
frequencies  using  very  small  data  sets  md  its  ability  to  produce  an  estimate  even 
when  the  correlation  matrix  is  not  of  full  rank  .  Additionally,  its  performance  in  a 
high  noise  environment  wjts  noteworthy. 

The  algorithm  also  proved  to  be  a  viable  method  for  data-adaptive  image 
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restoration.  The  nature  of  this  application  differs  greatly  from  that  of  spectral  estima¬ 
tion  in  that  the  data  sets  are  much  larger  and  generally  less  correlated.  Encouraging 
results  were  obtained  for  the  data-adaptive  noise  canceler  and  line  enhancer  with  only 
minor  modifications  of  the  same  algorithm  that  was  used  for  the  spectrum  estimates. 

B.  RECOMMENDATIONS  FOR  FUTURE  STUDY 

The  primaury  objective  of  this  thesis  was  to  develop  the  data-adaptive  iterative 
Toeplitz  algorithm  amd  obtain  experimental  results  to  evaduate  its  potential  for  use  in 
applications  that  involve  2-D  AR  parameter  estimation.  There  are  mamy  aspects  re- 
gau'ding  this  work  that  require  more  in-depth  study  and  several  applications  that  may 
benefit  in  the  use  of  this  adgorithm.  In  particular,  more  detailed  amalysis  of  model 
orders  amd  their  relation  to  the  algorithm’s  ability  to  form  spectral  estimates  from 
very  smadl  input  arrays  needs  to  be  cauried  out.  The  use  of  singulau^  vadue  decompo¬ 
sition  (SVD)  and  eigenvalue  amadysis  of  the  observed  data  may  provide  some  insights 
regarding  the  performauace  of  the  iterative  adgorithm.  Also,  further  investigation  as  to 
why  one  region  of  support  (ROS)  provides  greater  resolution  than  another,  may  yield 
information  that  would  permit  the  choice  of  am  optimum  ROS  prior  to  processing 
the  data.  Further  amadysis  for  choosing  vadues  of  the  forgetting  factor  A  amd  how  it 
relates  to  the  statistical  natme  of  the  data,  is  desirable  ais  well.  This  may  lead  to  a 
scheme  for  on-line  adjustment  of  A,  which  would  be  pairticulaxly  advantageous  when 
processing  images.  Other  future  work  could  explore  the  use  of  different  data  win¬ 
dowing  schemes  and  data-ordering  schemes  for  patssing  the  filter  mask  over  the  data. 
One  possible  method  that  has  not  been  investigated  is  an  expanding  square  starting 
at  one  comer  of  the  auray.  This  method  seems  well  suited  for  block-by-block  image 
processing.  Finally,  the  data-adaptive  iterative  Toeplitz  approximation  adgorithm  has 
the  potentiad  for  use  in  non-linear  applications,  such  as  the  identification  of  Volterra 
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